#' ---
#' title: "Robustness Check 3 (WVS7)"
#' author: "Gento Kato & Fan Lu"
#' date: "June 8, 2023"
#' ---
#' 
#' 
#' # Preparation 
#' 

## Clean Up Space
rm(list=ls())

## Set Working Directory (Automatically) ##
setwd(dirname(rstudioapi::getActiveDocumentContext()$path)); 

require(interflex)
require(estimatr)
require(texreg)
require(ggplot2)
require(ggrepel)
require(haven)

wvs7 <- readRDS("wvs7_japan_v1.rds")

cfmap <- list("(Intercept)"="(Intercept)",
              "edu2"="University Education",
              "cohortCohort II (18+ in 1975-1989)"="Cohort II (18+ in 1975-89)",
              "cohortCohort III (18+ in 1990-99)"="Cohort III (18+ in 1990-99)",
              "cohortCohort IV (18+ in 2000-09)"="Cohort IV (18+ in 2000-09)",
              "cohortCohort V (18+ in 2010-)"="Cohort V (18+ in 2010-)",
              "edu2:cohortCohort II (18+ in 1975-1989)"="University * Cohort II",
              "edu2:cohortCohort III (18+ in 1990-99)"="University * Cohort III",
              "edu2:cohortCohort IV (18+ in 2000-09)"="University * Cohort IV",
              "edu2:cohortCohort V (18+ in 2010-)"="University * Cohort V",
              "male"="Gender (Male)", 
              "cohortCohort II (18+ in 1975-1989):male"="Male * Cohort II",
              "cohortCohort III (18+ in 1990-99):male"="Male * Cohort III",
              "cohortCohort IV (18+ in 2000-09):male"="Male * Cohort IV",
              "cohortCohort V (18+ in 2010-):male"="Male * Cohort V",
              "incomecatMiddle" = "Income (Middle)",
              "cohortCohort II (18+ in 1975-1989):incomecatMiddle"="Income (Middle) * Cohort II",
              "cohortCohort III (18+ in 1990-99):incomecatMiddle"="Income (Middle) * Cohort III",
              "cohortCohort IV (18+ in 2000-09):incomecatMiddle"="Income (Middle) * Cohort IV",
              "cohortCohort V (18+ in 2010-):incomecatMiddle"="Income (Middle) * Cohort V",
              "incomecatHigh" = "Income (High)",
              "cohortCohort II (18+ in 1975-1989):incomecatHigh"="Income (High) * Cohort II",
              "cohortCohort III (18+ in 1990-99):incomecatHigh"="Income (High) * Cohort III",
              "cohortCohort IV (18+ in 2000-09):incomecatHigh"="Income (High) * Cohort IV",
              "cohortCohort V (18+ in 2010-):incomecatHigh"="Income (High) * Cohort V",
              "incomecatMissing" = "Income (Missing)",
              "cohortCohort II (18+ in 1975-1989):incomecatMissing"="Income (Missing) * Cohort II",
              "cohortCohort III (18+ in 1990-99):incomecatMissing"="Income (Missing) * Cohort III",
              "cohortCohort IV (18+ in 2000-09):incomecatMissing"="Income (Missing) * Cohort IV",
              "cohortCohort V (18+ in 2010-):incomecatMissing"="Income (Missing) * Cohort V",
              "workstatStudent/Housemaker/Part-Time" = 
                "Student/Housemaker/Part-Time",
              "cohortCohort II (18+ in 1975-1989):workstatStudent/Housemaker/Part-Time"=
                "Student/Housemaker/Part-Time * Cohort II",
              "cohortCohort III (18+ in 1990-99):workstatStudent/Housemaker/Part-Time"=
                "Student/Housemaker/Part-Time * Cohort III",
              "cohortCohort IV (18+ in 2000-09):workstatStudent/Housemaker/Part-Time"=
                "Student/Housemaker/Part-Time * Cohort IV",
              "cohortCohort V (18+ in 2010-):workstatStudent/Housemaker/Part-Time"=
                "Student/Housemaker/Part-Time * Cohort V",
              "workstatStudent/Part-Time" = 
                "Student/Part-Time",
              "cohortCohort II (18+ in 1975-1989):workstatStudent/Part-Time"=
                "Student/Part-Time * Cohort II",
              "cohortCohort III (18+ in 1990-99):workstatStudent/Part-Time"=
                "Student/Part-Time * Cohort III",
              "cohortCohort IV (18+ in 2000-09):workstatStudent/Part-Time"=
                "Student/Part-Time * Cohort IV",
              "cohortCohort V (18+ in 2010-):workstatStudent/Part-Time"=
                "Student/Part-Time * Cohort V",
              "workstatSelf-Employed/Full-Time/Managerial" = 
                "Self-Employed/Full-Time",
              "cohortCohort II (18+ in 1975-1989):workstatSelf-Employed/Full-Time/Managerial"=
                "Self-Employed/Full-Time * Cohort II",
              "cohortCohort III (18+ in 1990-99):workstatSelf-Employed/Full-Time/Managerial"=
                "Self-Employed/Full-Time * Cohort III",
              "cohortCohort IV (18+ in 2000-09):workstatSelf-Employed/Full-Time/Managerial"=
                "Self-Employed/Full-Time * Cohort IV",
              "cohortCohort V (18+ in 2010-):workstatSelf-Employed/Full-Time/Managerial"=
                "Self-Employed/Full-Time * Cohort V",
              "married"="Currently Married",
              "cohortCohort II (18+ in 1975-1989):married"="Married * Cohort II",
              "cohortCohort III (18+ in 1990-99):married"="Married * Cohort III",
              "cohortCohort IV (18+ in 2000-09):married"="Married * Cohort IV",
              "cohortCohort V (18+ in 2010-):married"="Married * Cohort V",
              "urban"="Urbanness (Current Residence)",
              "cohortCohort II (18+ in 1975-1989):urban"="Urbanness * Cohort II",
              "cohortCohort III (18+ in 1990-99):urban"="Urbanness * Cohort III",
              "cohortCohort IV (18+ in 2000-09):urban"="Urbanness * Cohort IV",
              "cohortCohort V (18+ in 2010-):urban"="Urbanness * Cohort V")

#'
#' # Analysis (Figure A4, Table A5)
#'

#+
##########################
## Analysis of 2019 WVS ##
##########################

## Models with Categorical Interaction
mA_wvs7 <- lm_robust(immigincrease~edu2*cohort+
                           male*cohort+
                           incomecat*cohort+
                           workstat*cohort+
                           married*cohort+
                           urban*cohort, 
                         data=wvs7, se_type="stata")
screenreg(list(mA_wvs7), 
          digits=3, include.ci=FALSE, single.row=TRUE,
          custom.model.names = c("Admit Foreign Workers"),
          custom.coef.map = cfmap)
texreg(list(mA_wvs7), 
       digits=3, include.ci=FALSE, single.row=TRUE,
       custom.model.names = c("Admit Foreign Workers"),
       custom.coef.map = cfmap, 
       custom.note = "%stars. Robust standard errors in parentheses.",
       booktabs = TRUE, dcolumn = TRUE, use.packages = FALSE, table=FALSE,
       file = "mtab_wvs7.tex")

mA_edu2_wvs7 <- 
  as.data.frame(  rbind(
    summary(lm_robust(immigincrease~edu2*relevel(cohort,levels(cohort)[1])+
                        male*relevel(cohort,levels(cohort)[1])+
                        incomecat*relevel(cohort,levels(cohort)[1])+
                        workstat*relevel(cohort,levels(cohort)[1])+
                        married*relevel(cohort,levels(cohort)[1])+
                        urban*relevel(cohort,levels(cohort)[1]), 
                      data=wvs7, se_type="stata"))$coefficients[2,c(1,4,5,6)],
    summary(lm_robust(immigincrease~edu2*relevel(cohort,levels(cohort)[2])+
                        male*relevel(cohort,levels(cohort)[2])+
                        incomecat*relevel(cohort,levels(cohort)[2])+
                        workstat*relevel(cohort,levels(cohort)[2])+
                        married*relevel(cohort,levels(cohort)[2])+
                        urban*relevel(cohort,levels(cohort)[2]), 
                      data=wvs7, se_type="stata"))$coefficients[2,c(1,4,5,6)],
    summary(lm_robust(immigincrease~edu2*relevel(cohort,levels(cohort)[3])+
                        male*relevel(cohort,levels(cohort)[3])+
                        incomecat*relevel(cohort,levels(cohort)[3])+
                        workstat*relevel(cohort,levels(cohort)[3])+
                        married*relevel(cohort,levels(cohort)[3])+
                        urban*relevel(cohort,levels(cohort)[3]), 
                      data=wvs7, se_type="stata"))$coefficients[2,c(1,4,5,6)],
    summary(lm_robust(immigincrease~edu2*relevel(cohort,levels(cohort)[4])+
                        male*relevel(cohort,levels(cohort)[4])+
                        incomecat*relevel(cohort,levels(cohort)[4])+
                        workstat*relevel(cohort,levels(cohort)[4])+
                        married*relevel(cohort,levels(cohort)[4])+
                        urban*relevel(cohort,levels(cohort)[4]), 
                      data=wvs7, se_type="stata"))$coefficients[2,c(1,4,5,6)],
    summary(lm_robust(immigincrease~edu2*relevel(cohort,levels(cohort)[5])+
                        male*relevel(cohort,levels(cohort)[5])+
                        incomecat*relevel(cohort,levels(cohort)[5])+
                        workstat*relevel(cohort,levels(cohort)[5])+
                        married*relevel(cohort,levels(cohort)[5])+
                        urban*relevel(cohort,levels(cohort)[5]), 
                      data=wvs7, se_type="stata"))$coefficients[2,c(1,4,5,6)]
  )  )
colnames(mA_edu2_wvs7) <- c("est","p","lci","uci")
mA_edu2_wvs7$dv <- "Admit Foreign Workers"
unique(sort(wvs7$univyr))
median(1943:1989); median(1975:1989); median(1990:1999); median(2000:2009); median(2010:2022)
mA_edu2_wvs7$x <- c(1966,1982,1994.5,2004.5,2016)

m_edu2_wvs7 <- mA_edu2_wvs7
m_edu2_wvs7$dv <- factor(m_edu2_wvs7$dv,levels=unique(m_edu2_wvs7$dv))
m_edu2_wvs7$lb <- factor(c("I","II","III","IV","V"), levels=c("I","II","III","IV","V"))
m_edu2_wvs7

## Interflex (without background variables)
set.seed(12345)
mxA_wvs7 <- interflex(Y = "immigincrease", D = "edu2", X = "univyr",
                          Z = c("male", "incomecat", "workstat", "married", "urban"),
                          data = wvs7, na.rm=TRUE, estimator = "kernel",
                          full.moderate = TRUE, bw=5, #CI=FALSE,
                          nboots = 1000, parallel = TRUE, cores = 3, theme.bw=TRUE)
# mxA_wvs7$figure

# save(mxA_wvs7, file="mx_wvs7.RData")
# load("mx_wvs7.RData")

tmpA <- as.data.frame(mxA_wvs7$est.kernel$`1`)
tmpA$dv <- "Admit Foreign Workers"
out <- tmpA
out$dv <- factor(out$dv,levels=unique(out$dv))

tmpA <- data.frame(x=mxA_wvs7$hist.out$mids,
                   y=c(mxA_wvs7$count.tr$`0`,mxA_wvs7$count.tr$`1`),
                   tr=rep(c("0","1"), each=length(mxA_wvs7$hist.out$mids)))
tmpA$dv <- "Admit Foreign Workers"
outh <- tmpA
outh$dv <- factor(outh$dv,levels=unique(outh$dv))
outh$tr <- factor(outh$tr, levels=c("1","0"))

#+ fig.width=5, fig.height=4
## Figure
movement7 <- 0.2
ggplot(out) + 
  geom_vline(aes(xintercept=1975), linetype=1, color="gray80") + 
  geom_vline(aes(xintercept=1990), linetype=1, color="gray80") + 
  geom_vline(aes(xintercept=2000), linetype=1, color="gray80") + 
  geom_vline(aes(xintercept=2010), linetype=1, color="gray80") + 
  geom_col(data=outh, aes(x=x,y=y*0.002,fill=tr), width=0.8) + 
  geom_hline(aes(yintercept=0+movement7), linetype=2, color="gray50") + 
  geom_ribbon(aes(x=X,ymin=`lower CI(95%)`+movement7,
                  ymax=`upper CI(95%)`+movement7),alpha=0.3) +
  geom_line(aes(x=X,y=TE+movement7)) +
  # geom_line(aes(x=X, y=`lower CI(95%)`+movement7),alpha=0.7, linetype=2) +
  # geom_line(aes(x=X, y=`upper CI(95%)`+movement7),alpha=0.7, linetype=2) +
  geom_errorbar(data=m_edu2_wvs7, aes(x=x, ymin=lci+movement7,ymax=uci+movement7), alpha=0.8, width=2) + 
  geom_point(data=m_edu2_wvs7, aes(x=x, y=est+movement7), size=1.5, color="gray10", alpha=0.8) + 
  geom_text(data=m_edu2_wvs7, aes(x=x, y=-0.23+movement7, label=lb), family="serif") + 
  facet_grid(.~dv) + 
  scale_fill_manual(values=c("black","gray70")) + 
  scale_x_continuous(breaks=c(1975,1990,2000,2010)) +
  scale_y_continuous(labels = function(y) round(y - movement7,2)) +
  coord_cartesian(ylim=c(-0.02,0.6)) + 
  labs(x = "Year of Turning 19 (Typical Year of Entering University)",
       y = "Conditional Coefficient with 95% Confidence Interval",
       caption = paste0("Histogram indicates the frequency of cases (black = university educated, gray = not).")) +
  theme_classic(base_family="serif") + 
  theme(plot.title = element_text(hjust=0.5),
        panel.background = element_rect(color="black", size=1),
        strip.text = element_text(size=11, face="bold"),
        legend.position="none")
ggsave("mx_cplot_wvs7.pdf", width=5, height=4)